'''
本示例用sympy，numpy，scipy库求解齐次线性方程组和非齐次线性方程组。
对于齐次线性方程组Ax = 0，用scipy库求解零空间的函数nullspace即可
对于非齐次线性方程组Ax = b，解的情况根据R(A)和R(A,b)的关系分为三种情况：
(1)R(A) = R(A,b) = n，有唯一解
(2)R(A) < R(A,b)，无解
(3)R(A) = R(A,b) < n，有无穷多解
x = pinv(A) @ b总是给出唯一数值解(x = A^(-1)b)，若是情况(2)，则给出的是最小二乘解,即Ax-b的模最小时的解；若是情况(3)，则给出最小范数解

本示例内容：
(1)求齐次方程组  x1 - 5x2 + 2x3 - 3x4 = 0; 的基础解系。其中向量的个数为n-r(A), n为变量个数，r(A)为矩阵A的秩
                5x1 + 3x2 + 6x3 - x4 = 0; 
                2x1 + 4x2 + 2x3 + x4 = 0
                
(2)求非齐次方程组  x1 + x2 - 3x3 - x4 = 1; 的解 （对应情况（3）,应获得最小范数解
                  3x1 - x2 - 3x3 + 4x4 = 4
                  x1 + 5x2 - 9x3 - 8x4 = 0

(3)利用numpy的lstsq()函数求超定方程组(方程数小于变量数) 0m + 1c = -1;的最小二乘解 ，相当于给定了x的观测值y，拟合经验函数直线，两个变量是斜率和截距
                                                    1m + 1c = 1;   
                                                    2m + 1c = 0.9;
                                                    3m + 1c = 2.1;
                                
'''

import numpy as np
from scipy.linalg import null_space, pinv
import sympy as sp
from matplotlib import pyplot as plt
# 定义齐次线性方程组的系数矩阵A
A = np.array([[1, -5, 2, -3],[5, 3, 6, -1],[2, 4, 2, 1]])
B = sp.Matrix(A)  # 将numpy数组转换为sympy矩阵
# 求解零空间，得到基础解系
print(f"矩阵A的秩为: {np.linalg.matrix_rank(A)}") 
print(f"采用sympy库求解零空间的结果为: {B.nullspace()}")  # sympy的nullspace方法可以直接求解零空间
print(f"采用scipy库求解零空间的结果为: {null_space(A)}")  # scipy的null_space方法可以直接求解零空间
#两个显示结果不同，但是等价的；scipy库返回的是标准正交基，sympy库返回的是分数形式

# 定义非齐次线性方程组的系数矩阵A和常数向量b
A = np.array([[1, 1, -3, -1],[3, -1, -3, 4],[1, 5, -9, -8]])
b = np.array([1, 4, 0]).reshape(3,1)
#定义增广矩阵(A,b)
A_augmented = np.hstack([A,b])  # 将A和b水平堆叠(在axis=1上操作)形成增广矩阵
print(f"A的秩为{np.linalg.matrix_rank(A)}，增广矩阵的秩为{np.linalg.matrix_rank(A_augmented)},变量个数为{A.shape[1]}")
x = pinv(A) @ b  # 求解非齐次线性方程组的最小范数解
print(f"最小范数解为{x}")

x = np.array([[0,1],[1,1],[2,1],[3,1]]).reshape(4,2)
y = np.array([-1,0.2,0.9,2.1]).reshape(4,1)  
m,c = np.linalg.lstsq(x,y, rcond=None)[0]
print(f"最小二乘解为{np.linalg.lstsq(x,y, rcond=None)[0]}")  # 求解超定方程组的最小二乘解
plt.rcParams['font.family'] = 'SimHei'  # 设置为黑体（或其他支持中文的字体）
plt.rcParams['axes.unicode_minus'] = False  # 正确显示负号
plt.figure(0)
plt.plot(x[:,0],y.reshape(-1),'o') # 绘制散点用'o' marker
plt.plot(x[:,0],m*x[:,0]+c,'r')
plt.legend('原始数据', '拟合直线')
plt.show()
